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I.  INTRODUCTION 


The  objective  of  this  work  is  to  develop  an  accurate  and  efficient 
numerical  procedure  for  solving  the  unsteady  Navier-Stokes  equations  tc. 
describe  transient  spin-up  flow  occurring  in  a  cylindrical  container  when 
it  is  suddenly  rotated  about  its  longitudinal  axis.  Knowledge  of  this 
internal  flow  is  needed  to  design  gun-launched  projectilos  which  carry 
smoke/ incendiary  agents  or  chemical  payloads.  Liquid  payloads  enhance 
spin  decay  of  project i lcs1 and  their  presence  can  produce  flight  dynamic 
instabilities  as  a  result  of  resonance  between  the  projectile  nutational 
motion  and  inertial  oscillations  in  the  rotating  liquid3.  Prom  a 
computational  viewpoint  this  problem  is  instructive  because  it  is  an 
example  of  a  class  of  internal  flow  problems  for  which  computational 
experiments  can  uncover  details  of  the  flow  that  cannot  be  easily 
visualized  or  measured  experimentally. 

The  results  presented  here  demonstrate  that  a  predictor-corrector 
mult iplc- iteration  (PCM!)  technique  developed  by  Rubin  and  LinM  for 
solving  steady  three-dimensional  boundary  region  problems  can  be  success¬ 
fully  adapted  to  solve  the  unsteady  Navier-Stokes  equations.  In  the 
present  approach  this  racthoJ~Ts  combined  with  ie  Gauss-Seidel  procedure ’’ 
and  grid  stretching  transformations  to  produce  an  accurate  and  efficient 
numerical  procedure  for  describing  the  spin-up  process.  Calculations 
with  the  PCMI  method  have  been  performed  for  spin-up  from  rest  and  spin- 
up  from  an  initial  state  of  solid-body  rotation;  in  both  types  of 
problems  inertial  oscillations  have  developed  in  the  rotating  liquids. 
Numerical  results  have  been  obtained  for  a  range  of  cylinder  aspect 
ratios,  n,  from  0.3  to  4.4  and  a  range  of  Reynolds  numbers  from  2 1 S  to 
50,000.  Calculations  performed  for  five  test  problems  are  consistent 


lF.  H.  Wedemt'yer,  "The  Unsteady  Flout  Within  a  Spinning  Cylinder," 

.7,  Fluid  /■k'oh . ,  Vol.  20,  Ft.  d,  1064,  pp.  383-399;  aUo  see  BRl  Report 
1262,  Aberdeen  Moving  Ground,  ND,  AD  431646,  Oot.  1963. 

*V.  W.  Kitchens,  Jr.,  N.  Gerber  and  R.  Sedney,  " Spin  Decay  of  Liquid- 
Filled  Project! lee, "  J.  Si^ace  craft  and  Rockets,  Vol.  IS,  No.  6,  Nov- Dec 
1976,  pp.  348-3S4. 

3A‘.  Stewartson,  "On  the  Stability  of  a  Spinning  Top  Containing  Liquid 
J,  Fluid  Me ah,  Vol.  S,  Pt.  4,  Sept.  19S9,  pp.  S77-S92. 

US.  G.  Rubin  and  T.  C.  Lin,  "A  Nienerical  Method  for  Three-Dimensional 
Viscous  Floui:  Application  to  the  Hypereonio  Leading  Edget " 

J.  Comp.  Fhye.t  Vol .  9,  1972 ,  pp.  339-364, 


5M.  G.  Salvadori  and  M.  L.  Baron,  Nxsmrioal  Methods  in 
Prentice-Hall,  Inc.,  Engleixtod  Cliffs,  N.  J,,  1961. 
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with  previous  computations6* 8  and  experimental  measurements6 •  '*, 
Numerical  results  huve  also  been  used  to  quantify  the  flow  in  the  likmun 
(or  endwall)  boundary  layers  during  spin-up  and  thus  develop  an 
appropriate  "compatibility  condition"  for  use  in  Wedemeyer's  model1  for 
spin-up  from  rest.  Although  the  results  are  not  discussed  here,  the 
PCMI  procedure  has  been  used  to  obtain  spin-up  flow  in  the  annulus 
between  finite-length,  concentric  cylinders.  Neiticl10  also  used  the 
PCMl  computer  program  to  study  the  onset  and  tenoral  development  of 
fluid  dynamic  instabilities  during  spin  down  in  a  cylinder. 


II.  GOVERNING  EQUATIONS  AND  BOUNDARY  CONDITIONS 

The  calculations  employ  a  finite-difference  analog  of  the  unsteady 
axisymmetric  Nav ier-Stokes  equations  formulated  in  cylindrical  coordinates 
ir.o.z).  The  equations  are  expressed  in  terms  of  stream  function,  f, 
vorticity,  C.  and  circulation,  y,  instead  of  velocity  and  pressure  in 
order  to  simplify  the  numerical  procedure.  In  dimensionless  variables 
the  governing  equations  are 


v2r  -  vr/r  ■  rc,  (1) 

t.  ♦  uc  ♦  wc  -  u c/r  -  Syy./r1  ■  (l/Re)|V2c  ♦  C_/r  -  c/r2],  (2) 

i  r  z  z  r 

Yt  *  UYr  *  ”  Yr/rl;  (3) 

where  the  subscripts  denote  partial  differentiation  and 

Re  «  Oa2/v,  (4) 


.  bhmi-Varnas ,  V.  V.  Foul  in ,  S.  Piacsek  and  5.  M.  Lee,  "Numerical 
'lutiow  ami  Laser  !\?pp U} r  Measurements  of  Spin- Up,  "  J,  Fluid  Mech. , 
i  l.  86,  Pt .  4,  1978,  pp.  809-639. 

7W.  R.  Briley,  "Time  Dependent  Notating  Flou  in  a  Cylindrical  Container," 
PhD  Dissertation,  The  University  of  Texas  at  Austin,  1968,  University 
Microfilms,  Inc,,  69-617.1. 

8V.  R.  Briley  <vid  H.  A.  Valle,  "A  Nienerical  Study  of  Time-Dependent 
Rotating  Floo  in  a  Cylindrical  Container  at  Lent  ami  Moderate  Reynolds 
Nxerters," Proc.  Pnd  Inti.  Conf.  on  Num.  Heth.  Fid.  Dyn. ,  Lecture  Notes 
in  Physics,  Vol.  6,  Springer- Verlag,  1970,  pp.  377-384. 

9V.  B.  Uatkine  and  R.  G.  Hussey,  " Spin-Up  From  Rest  in  a  Cylinder , " 

Phys.  of  Fluids,  Vol.  20,  No.  10,  Pt.  1,  1977,  pp.  1696-1604. 

1  °J.  P.  Neitzel,  Jr.,  " Centrifugal  Instability  of  Decelerating  Svirl 
Flou  urithin  Finite  and  Infinite  Circular  Cylinders,"  PhD  Dissertation, 
The  John  Hopkins  University,  Baltimore,  HD.  1979. 


V‘ 


Y  c  rv , 


(5) 


(<*) 


|  with  the  axisymmctric  stream  function  defined  so  that 

i 

;■  u  «  t  /r  and  w  •  -  Y  /r. 

4  s  r 


(7) 


(8) 


t  The  hkman  number  based  on  ha  If -height  is  related  to  Re  by 

t 

!  l.k  »  v/(llc?)  =  l/(a2Re). 


k 

4 

*. 


The  stream  function-vorticity-circulation  formulation  yields  an  elliptic 
I’D!;,  liquation  (1),  and  two  parabolic  1'Dlis,  liquations  (2)  and  (5),  which 
are  coupled.  The  boundary  conditions  impose  additional  coupling  between 
y  and  r. . 

The  nondimensional  variables  used  here  arc  formed  by 
r  =  R/a,  2  »  Z/a,  t  *  UT, 

u  «  U/(P.a),  v  -  V/  (wa) ,  w  »  K/(na),  (9) 

t  «  «/(}?»-*),  y  *  I7(U  a  2 ) ,  r.  »  Z/fl. 

The  initial  conditions  for  spin-up  are 

Y  3  e  0 ,  >  =  u.r2/0  for  t  <  0,  (10) 

where  is  the  initial  cylinder  rotation  rate.  Computational  efficiency 
and  resolution  arc  improved  by  employing  a  symmetry  boundary  conditior 
at  the  cylinder  mid-plane,  r  *  «.  This  effectively  halves  the  number 
of  grid  points  required.  The  boundary  conditions  for  t  >  0  are 

Y(t.  0,  2)  «  y ( t ,  0,  2)  ■  t(t,  o,  I)  -  0,  (11a) 

nt,  1,  2)  *  0,  Y(t,  l,  2)  «  1,  c  ( t  ,  1,  Z)  -  Yrr  (t,  1,  2).  (lib) 
t(t,  r,  0)  «  0,  y (t ,  r,  0)  ■  r2,  t(t,  r,  0)  ■  r»  0)/*.  (He) 

f(t,  r,  a)  ■  ;(t,  r,  a)  «  0,  Yj(t,  r,  a)  «  0. 


(lid) 


The  boundary  conditions  tor  vortlolty  along  t ho  sidewall  and  i-ndwull, 
Equations  (lib)  and  (IK),  are  derived  from  Equations  (7)  and  (8)  by 
imposing  the  no-slip  conditions  for  velocity.  Figure  1  illustrates  the 
coordinate  system  and  boundaries  used  in  the  numerk  1  calculations 
for  spin -up. 

During  the  spin-up  process  there  are  viscous  regions  near  the  side- 
wall  and  endwalls  which  become  very  thin  as  Re  becomes  larger  the  1000 
or  so,  necessitating  a  fine  grid  to  resolve  the  boundary- layer  type 
phenomena  along  these  walls.  Analytical  coordinate  transformations  are 
used  to  optimize  the  grid  point  distribution  and  transform  a  nonuniform 
grid  in  the  physical  plane  into  an  equal  I y- spaced  grid  in  the  computa¬ 
tional  plane.  Transformations  based  on  the  work  of  Roberts’1, 


8  »  In  ((&  ♦  r)/(b  -  r) 1/ln  ( (S  ♦  l)/(b  -  1)),  (12a) 

n  ■  I  ♦  In  |  (e  ♦  2/u  -  l)/(c  -  i/a  ♦  l))/ln  ( (c  ♦  l)/(c  -  1)|,  (12b) 

-1/2 

with  b  »  (1  *  d)  *"  and  c  »  (1  -  e)  are  selected.  These  trans- 
format  ■  mis  are  particularly  suited  for  problems  where  thin  viscous 
regions  lie  along  one  boundary  in  each  of  the  coordinate  directions. 

Values  of  d  and  e  (0  <  d  <  1,  0  *  e  <  1)  are  specified  to  group  a  large 

fraction  of  the  grid  points  (typically  1/2  or  so)  into  the  sidewall  and 

endwall  viscous  regions  where  large  velocity  gradients  arc  present;  as 
d  and  e  *  0  the  grid  point  spacing  becomes  finer  near  the  sidewall  and 
endwall,  respectively.  Figure  2  shows  a  typical  nonuniform  grid  point 
distribution  in  the  physical  plane  produced  with  these  transformations,  sec 
Equation  (12).  using  equally-spaced  grid  points  in  the  computational 
plane.  The  complete  set  of  transformed  equations  and  boundary  conditions 
arc  giver,  in  Appendix  A. 


III.  NUMI.RIC.Al.  PROCF.lXJRi. 

Many  methods  have  been  used  by  previous  investigators  to  solve  the 
stream  function-vorticity  fonn  of  *hc  Navier-Stokes  equations.  Perhaps 
the  most  popular  technique  i  ..  to  combine  the  alternating-direction 
implicit  (ADD  method’2  for  the  t  and  y-equetions  with  cither  an  ADI  or 


1  0 .  Roberto,  "Comp.utational  Meoheo  for  Boundary- Layer  Problem o," 

Proa.  Pnd  Int  i.  Conf.  on  than.  Hath.  Fid.  Dyn. ,  Lecture  Noteo  in  Phyeico, 
Vol.  8,  Springer  Vorlag,  1970,  pp.  171-177. 

32P.  V.  Pcaaman  and  H.  H.  Raohford,  Jr.t  nThe  Nunerioal  Solution  of 
Parabolic  and  Elliptic  Differential  Eqwitione , "  J.  SIAN ,  Vol.  3,  No. 

1 ,  195t>,  pp.  P.8-41. 
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successive  over-relaxation  (SOR)  method13  for  the  t-equatlon.  This 
report  describes  an  efficient  alternative  procedure  for  solving  the  t, 
and  * -equations,  namely,  the  semi-implicit  PCM1  method.  To  the  best  of 
the  author's  knowledge  this  is  the  first  time  this  method  has  been  used 
to  soivo  the  Navier-Stokes  equations.  It  represents  a  compromise  in 
approach  between  the  implicit  ADI  scheme  and  explicit  schemes  used  by 
other  authors.  In  the  nresent  application  the  V-equation  is  solved  by 
the  Gauss-Scidel  method’;  the  SOR  method  was  used  for  test  calculations, 
but  it  did  not  speed  up  the  overall  procedure. 

In  the  PCMI  method  the  solution  is  advanced  to  a  new  time  level  in  a 
single  time  step  At  as  opposed  to  the  two  half-time  steps  required  in 
one  cycle  of  the  ADI  method.  It  is  implicit  in  the  radial  direction; 
the  solution  requires  only  the  inversion  of  a  tridiagonal  matrix  for  each 
row  in  the  computational  grid.  The  ADI  procedure,  on  the  other  hand, 
requires  that  a  tridiagonal  matrix  be  solved  for  each  row  in  the  first 
half-time  step,  and  then  for  each  column  in  the  second  half-time  step, 
leading  to  a  longer  computation  time  per  full  time  step. 


A  symmetry  boundary  condition,  such  a>  yt  3  0  along  the  cylinder 
mid-plane,  is  easy  to  implement  in  this  method  since  all  flow  gradients 
in  the  :-direction  are  approximated  by  prediction  and  subsequent 
correction  in  this  time-iteration  technique.  This  approach  eliminates 
the  cross  coupling  of  grid  points,  thus  reducing  the  size  of  the  inver¬ 
sion  matrices  and  decreasing  the  computer  time  required.  The  iteration 
procedure  allows  the  boundary  vorticity  to  converge  and  also  allows  the 
nonlinear  terms  to  be  approximated  and  then  corrected,  giving  a  more 
accurate  simulation  of  the  nonlinear  coupling  between  equations. 

Central  difference  formulae  are  used  for  all  spatial  derivatives  at 
interior  points,  avoiding  false-diffuslon  effects  introduced  by  upwind 
difference  schemes.  Temporal  derivatives  are  approximated  by  second- 
order  accurate  one-sided  difference  formulae  involving  three  time  levels. 
The  following  finite-difference  representations  are  used  for  the  t- 
equation: 


(1/2A6) 


[>♦»  -  1 
L  .  j*l.  k  *♦!.  j-l.k J  • 


*66 


(1/AflO  t 


[ 


m*l 

i*l ,  j ♦ 1  .  k 


-  2; 


m*l 

1*1,  j.k 


ct  -  (i/2At)  ['c"*! j  k  -  4<4j,k  4  ci-i,j,k] 


i-i.k  J  • 

(ISa) 

Ci*l.  j-l.kj 

(13b) 

Ci-l.j.k]  * 

(13c) 

1 *D.  Young,  "Iterative  Method*  for  Solving  Partial  Different*  Equation* 
of  Elliptic  Type,"  Trane.  Amr.  Math ,  Soo . ,  Vol.  7E ,  1964 ,  pp.  92-111. 
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S  “  (1/2An)  [ci*l.j.k*l  *  CiM,j,k-i] 


tn-  (l/2An)[t"MJ(ko  -  ^uuk.x]  . 


(13d) 
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whore  superscript  n  denotes  the  tine- iteration  number.  The  expressions 
for  the  v-equation  are  identical  in  fona  to  Equations  (15).  The  finite- 
difference  representations  tor  the  f -equation  are: 

Ye  -  U/2A0)[t^|JM  k  -  Tjllj-Uk]’ 

V  "  (l/Ar)  [Ti*t,j*l,k  '  2li*l,;},k  *  tis!.j-l,k]  •  (l4b) 

Yn  “  (l/2An)  [Vl,j,k*l  ‘  Ti*l,j,k-l]  *  (l4c) 

Tnn  "  (1/An,)  [*"♦!, j,k»I  *  2Ti*I,j,k  *  Ti*l,j,k-l]  ;  (14d) 

where  superscript  n  denotes  the  t-iteration  nunber.  Figure  5  shows  a 
finite-difference  stencil  for  this  scheme  with  all  derivatives  evaluated 
at  (i*l,j,k).  The  above  finite-difference  approximations  arc  compact, 
involving  only  one  time  level.  This  leads  to  prograaming  simplifications 
and  shorter  run  tines.  The  complete  set  of  finite-difference  equations 
is  given  in  Appendix  B.  These  finite-difference  expressions  insure 
truncation  errors  for  interior  points  of  0(At2,  AO2,  An2).  The  overall 


.rwnif. 


accuracy  of  he  method  depends  on  the  treatment  of  the  boundary  con¬ 
ditions  and  this  will  be  discussed  later. 

The  numerical  procedure  is  easy  to  implement.  For  the  predictor 
step,  or  first  iteration  of  each  time-iterative  cycle,  terms  in  the 
difference  equations  with  superscript  0  are  approximated  by  a  Taylor 
series  to  0(At  ’) : 


F° 

»♦! ,  j  ,k 


31- . 


i  •  3  • 


(15) 


During  the  first  two  time*  steps  extrapolations  of  O(At)  and  0(At?)  are 
used,  respectively.  The  use  of  Equation  (IS)  reduces  the  number  of 
iterations  required  to  achieve  accuracy  and  stability;  these  advantages 
must  be  weighed  against  possible  storage  problems  caused  by  the  addi¬ 
tional  planes  of  data  needed  for  the  extrapolation.  The  influence  of 
the  extrapolation  procedure  on  stability  and  iteration  convergence  has 
been  discussed  by  Rubin  and  tin14. 

After  extrapolating  guesses  for  v®,  and  at  time  (i*T),  the 
wall  vorticity  is  determined  using  the  Y"-values  and  Equations  (lib)  and 
(He).  The  manner  in  which  this  is  carried  out  deserves  special  comment, 
since  it  can  often  have  a  strong  influence  on  iteration  convergence. 

In  the  present  calculations  we  adopt  a  first-order  form  for  the  wall 
vorticity  boundary  conditions,  expressed  in  transformed  coordinates  as 


cj:|  -  2(*r)2T^,/tA6;)  ♦  O(AB). 


(lba) 


m*  1 

C  * 
n»o 


2(n.)?  y*+|/(rAn2)  ♦  O(An); 


(10b) 


where  subscript  w*l  represents  the  grid  point  adjacent  to  each  respec¬ 
tive  wall  point.  This  first-order  forrm  is  used  since  it  is  known  to  have 
the  least  adverse  effect  on  iteration  convergence  and  it  often  gives 
results  essentially  identical  to  higher-order  accurate  expressions1 . 

The  numerical  procedure  used  here  appears  to  he  fully  compatible  with 
second-order  accurate  expressions  for  wall  vorticity,  based  on  the 
results  of  test  calculations  for  «  »  1,  Re  ■  1000.  The  boundary  values 
for  >  along  the  midplane  are  determined  from  the  one-sided  difference 
expression 


14C,  V.  Kitchens,  Jr.,  "Separation  and  Reattaohmnt  Near  Square 
Protuberances  in  Lou  Reynolds  Nverber  Couette  Flou,  "  BRL  Report  169b, 
Aberdeen  Proving  Ground,  MD,  AD  773663,  Jan.  1974. 
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where  M-l  and  M-2  represent  the  first  and  second  points,,  respectively, 
adjacent  to  n«l. 

With  all  boundary  values  now  approximated,  the  difference  equations 
for  and  tm*  are  solved  with  the  PCMl  method  using  the  m-iteratc 
values  to  form  the  coefficients  of  the  nonlinear  terms.  The  calculations 
start  along  the  row  of  points,  M-l,  adjacent  to  the  raidplunc  and  work 
downward  toward  the  endwall.  The  derivatives  in  the  ti-di  recti  on  arc 
treated  implicitly,  thus  requiring  the  solution  of  a  tridiagonal  system 
of  equations  along  each  successive  row.  T^e  (a*l)-iterate  values  at 
li*l,j,Ml)  are  used  to  approximate  derivatives  in  the  n-dircction  at 
(i+l.j.k)  as  soon  as  they  become  available;  see  liquations  (13d)  and 
llSc). 

At  the  end  of  each  iteration  cycle  for  c  and  y,  the  difference  form 
of  the  stream  function  equation  is  solved  iteratively  by  the  Gauss-Seidel 
technique.  The  solution  is  obtained  by  starting  at  the  interior  grid 
point  adjacent  to  6»n«0  and  sweeping  first  in  6  and  then  n,  making  use 
of  updated  values  as  soon  as  they  become  available;  sec  liquations  (14). 
Convergence  is  assumed  when 

Max  .  .1 Tn+1  -  Yni  <  c  .  (18a) 

J  »  *  1 

This  is  typically  achieved  in  2-3  iterations  with  Cj  *  lxlo'7.  The 
converged  values  for  V  arc  now  used  to  update  the  boundary  values  for  t 
and  repeat  the  iteration  process  for  y  and  c.  The  iteration  process  is 
assumed  to  converge  when  both 

MaXj  j,  I  y**1  -  y*\  <  c2  ,  (18b) 

Max^  U**1  -  Cml  <  .  (18c) 

It  typically  requires  2-3  iterations  to  satisfy  Equations  (18b,  c)  with 
c2  *  r3  *  1*1°  >  at  very  early  time  5-10  iterations  are  generally  needed 

due  to  the  severe  flow  unsteadiness  caused  by  the  impulsive  start  and 
the  subsequent  inaccuracy  of  the  extrapolated  guesses. 


IV.  STABILITY  PROPERTIUS  01  NUMERICAL  I'ROCIilHJRI: 


Rubin  and  Lin1'  have  analyzed  the  interior  point  stability  of  the 
PCM  l  method  for  a  model  linear  equation.  Their  analysis  shows  that  the 
method  is  conditionally  stable  and  that  with  one  or  more  iterations  the 
stability  criterion  is  independent  of  Re.  On  the  basis  of  their  results 
the  appropriate  stability  criterion  for  our  calculations  is 

At  s  Min.  .  I  An/I  n  /Re  ♦  B  n  tVrl  ] .  09) 

J  r  r  p 

The  term  in  the  denominator  containing  Re  results  from  the  coordinate 
transformation  for  z;  this  term  vanishes  if  an  equally-spaced  grid  is 
used  in  the  z-direction.  liquation  (ID)  must  be  applied  cautiously  since 
the  governing  equations  are  actually  non-linear  and  the  boundary  condi¬ 
tion  treatment  has  not  been  included  in  the  analysis. 

Numerical  tests  were  conducted  to  assess  the  effect  of  violating  the 
above  stability  criterion.  These  tests  were  conducted  for  a  *  1  with 
Re  *  1000,  974,!  and  SO, 000,  using  three  combinations  of  grid  sizes  for 
each  Re  and  several  different  values  of  the  transformation  parameters 
d  and  e,  The  results  show  that  numerical  stability  is  always  achieved 
when  Equation  (19)  is  satisfied.  In  some  cases  the  PCMI  calculations 
remain  stable  with  At  as  large  as  1 S0\  of  the  maximum  allowable  value. 

In  general,  the  calculations  show  that  satisfying  Equation  (19)  is 
sufficient,  hut  not  necessary,  for  numerical  stability.  Equation  (19) 
was  satisfied  at  each  time  step  in  the  illustrative  examples  to  he 
discussed  next. 


V.  COMPARISON  WITH  PREVIOUS  WORK 

The  present  method  has  been  used  to  treat  the  problems  of  spin-up  from 
rest  and  spin-up  from  an  initial  state  of  solid-body  rotation.  K'e 
compare  our  results  with  those  of  Warn-Varnas  et  al.(>  for  the  latter 
problem.  They  used  an  ADI  technique  coupled  with  a  scheme  developed  by 
Kill  lams'1'  to  solve  the  velocity-pressure  form  of  the  Navicr-Stokes 
equations.  In  their  calculations  they  differenced  the  governing  equations 
directly  on  a  stretched  grid  instead  of  transforming  to  new  coordinates 
as  is  done  here.  Their  computations  were  verified  by  measurements  taken 
with  a  laser  doppler  vclocimotcr  (U)V)  system. 

Figure  4  shows  a  comparison  of  the  present  calculations  with  results 
from  Reference  6  (their  Figure  13b).  The  comparisons  are  expressed  in 
terms  of  their  quantity  called  "zonal  velocity"  (ordinate  in  Figure  4) 


l5(7.P.  Uilliams,  "Numerical  Integration  of  the  Three-Dimensional  Navier- 
Stokee  Equations  for  Incompressible  Flw,”  J.  Fluid  Mech. ,  Vol.  37 , 
1969,  pp.  ?27-?b0. 
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Figure  4.  Inertial  osci nations  during  spin-up  from 
a  previous  state  of  rigid-body  rotation. 
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which  is  a  sealed  non -di mens ionu 1  angular  velocity.  The  results  are 
shown  at  r  «  0.25  on  the  cylinder  symmetry  plane  for  a  case  with  u  * 
0.5182,  Re  ■  7554  and  ih  ■  0.81820.  The  inertial  oscillations  excited 
by  the  sudden  increase  In  cylinder  rotation  rate  are  cloarly  predicted 
in  both  computations  and  are  in  fairly  good  agreement  with  experimental 
measurements.  Both  of  these  numerical  results  appear  to  be  within  the 
experimental  uncertainty  associated  with  these  data,  according  to  the 
error  analysis  presented  in  Reference  6.  Comparisons  for  several  other 
positions  in  the  cylinder  (not  shown  here)  give  similar  agreement  for  both 
the  decay  of  the  zonal  velocity  and  the  amplitudes  and  phases  of  the 
inertial  oscillations.  The  computation  time  end  number  of  grid  points 
used  to  obtain  the  numerical  results  in  Reference  6  are  not  stated.  The 
PCM l  method  required  51,5s  of  CPU  time  on  a  CDC  7600  computer  using  a 
stretched  (d  »  0.5,  e  ■  0.1)  41  x  21  (r-z)  grid  with  600  time  steps 
(At  ■  0.065).  Approximately  2-3  iterations  were  required  per  time  step. 

The  problem  of  spin-up  from  rest  has  been  emphasized  in  the  present 
work  because  of  its  application  to  liquid-filled  projectiles.  This 
problem  is.  by  its  nature,  nonlinear;  the  previous  problem  can  be 
linearized  for  small  0  -  0|.  Figure  5  compares  results  for  zonal  velocity 
for  spin-up  from  rest  with  those  for  spin-up  from  a  previous  state  of 
rigid-body  rotation.  The  comparison  is  made  at  r  »  0.90  for  two  values 
of  z,  illustrating  the  axial  structure  present  in  these  oscillations. 

For  •  0,  results  for  both  values  of  z  indicate  that  the  frequency  of 
the  dominant  inertial  mode  increases  with  time  at  the  early  times  shown 
in  Figure  5;  the  amplitude  of  these  oscillations  damps  rapidly  and  cannot 
be  detected  for  t  >  40. 

Comparisons  have  been  made  with  computations  of  Briley7  and  Briley 
and  Walls*'  for  spin-up  from  rest.  They  studied  this  problem  for  low  Re 
using  the  API  technique  to  solve  the  stream  funct ion-vort ic ity  form  of 
the  Navier-Stokes  equations.  Figure  6  compares  values  of  rotational 
volume  flow  rate, 

Q  -  ( l/o)  f1  f2"  v  dzdr,  (20) 

o  o 

for  two  cases.  The  quantity  Q  can  be  used  to  obtain  a  measure  of  the 
spin-up  time.  Briley  and  Walls  used  a  uniform  grid  that  became  re¬ 
strictive  at  moderate  Re  due  to  the  small  thickness  of  the  cndwall 
boundary  layers;  they  obtained  results  for  Re  as  large  as  1167.  Our 
calculations  appear  to  be  in  good  agreement  with  all  of  their  results 
for  spin-up.  NeitzePs  comparisons  for  spin-down10,  however,  showed 
only  qualitative  agreement  with  Briley  and  Walls'  results  for  Re  ■  1167. 
The  observed  differences  are  thought  to  be  due  to  grid  size  effects. 

The  present  computations  have  also  been  compared  with  LDV  measurements 
taken  by  Watkins  and  Hussey9.  Figure  7  presents  comparisons  of  azimuthal 
velocity  along  the  cylinder  mid-plane  at  four  instants  during  spin-up  for 
a  case  with  a  ■  1.S1S,  Re  «  5076.  Figure  8  shows  similar  comparisons  for 
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Figure  S.  Comparison  of  inertial  oscillations  during 
spin-up  from  rest  and  from  a  previous 
state  of  rigid-body  rotation. 
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I'iRUre  6.  Rotational  volume  flo«  rate  for  apin-up  from  test 
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1  ■  l.  Re  ■  9741.0.  The  size  of  the  symbols  used  to  plot  the  experi¬ 
mental  data  in  Figure*  7  and  8  annroximatel v  renresent*  the  size  of  the 
error  bars  that  should  be  attached  to  these  data.  The  calculations  in 
figure  8  used  a  .11  x  21  grid  with  d  ■  c  ■  0.10  and  required  2745  time 
steps  with  At  «  0.10  to  reach  t  ■  274.5.  Approximately  1-2  iterations 
were  required  per  time  step  to  satisfy  liquations  (18).  It  is  interesting 
to  note  that  the  PCM1  solution  required  only  69*  of  CPI)  time  on  the  CPC 
7t»00,  whereas  the  experimental  spin-up  process  depicted  in  figure  8 
required  ISOs. 

The  results  shown  in  figures  6,  7  and  R  are  representative  of  the 
"core"  flow  in  Wedemeyer's  model  of  spin-up  from  rest  and  they  can  be 
predicted  fairly  well  using  that  model;  the  accuracy  of  the  prediction 
increases  as  Re  increases.  However,  Wedcmcyer's  model  says  very  little 
about  the  flow  in  the  F.kman  layers,  in  the  corner  region,  and  along  the 
sidewall.  These  phenomena  will  be  discussed  next. 


VI.  TRANSIENT  PHENOMENA  AT  I  AK1.Y  TIME 

The  calculations  were  used  to  examine  the  details  of  the  spin-up 
flow  in  the  endwall  f.kman  layers  and  investigate  transient  reversed  fiow 
regions  ,wit  develop  and  then  decay  along  the  sidewall  during  the  first 
few  rotations  after  the  impulsive  start.  The  latter  phenomena  are 
illustrated  in  figures  9a-d.  Instantaneous  streamlines  arc  shown  for 
the  Watkins  and  Hussey  case  with  <>  *  1 ,  Re  »  9741.6,  based  on  calculations 
performed  with  a  41  x  41  grid  with  d  ■  e  ■  0.1  and  At  *  0.05.  The  calcu¬ 
lations  predict  the  development  of  several  weak  reversed  flow  regions 
in  the  meridional  flow  along  the  sidewall.  A  single  reversed  flow 
region  has  formed  near  the  corner  by  t  *  6  (see  figure  9b);  two  such 
regions  have  formed  along  the  sidewall  by  t  »  15  and  there  are  four 
present  by  t  •  20  (sec  figure  9c).  These  weak  reversed  flow  regions 
"collapse"  in  the  next  half-rotation  or  so  (sec  figure  9d  for  t  »  24) 
and  do  not  redevelop  for  t  s  24.  Grid  convergence  studies  for  this  case 
show  that  the  quantitative  results  in  figures  9  are  sensitive  to  grid 
si:c;  nevertheless,  the  qualitative  presence  of  the  reversed  flow  regions 
is  predicted  for  calculations  with  11  x  11,  21  x  21  and  41  x  41  grids. 

The  transient  reversed  flow  regions  do  not  develop  in  calculations 
for  (i  •  1,  Re  <  1000,  possibly  because  of  the  large  amount  of  viscous 
dissipation.  At  higher  Re  there  is  less  viscous  dissipation  present  and 
inertial  effects  become  more  pronounced.  At  very  early  time  the  inertial 
oscillations  are  confined  to  a  thin  layer  of  rotating  fluid  along  the 
sidewall.  Fluid  particles  near  the  endwall  are  accelerated  radially 
outward  in  a  spiral  motion  as  the  F.kman  layer  develops.  These  particles 
overshoot  their  "equilibrium  radial  position"  before  they  turn  upward 
from  the  edge  of  the  Ekman  layer  near  the  corner.  The  reversed  flow 
regions  that  dovelop  along  the  sidewall  are  apparently  linked  to  the 
inertial  oscillations  developed  as  swirling  fluid  particles  travel  up¬ 
ward  along  the  sidewall  and  begin  to  migrate  radially  inward.  As  Re 


increases,  the  calculations  predict  that  both  the  inertial  oscillutions 
und  reversed  flow  regions  become  more  pronounced.  Figure  10  shows 
Instantaneous  streamlines  for  o  »  1,  Re  »  50,000  and  t  «  15  ohtained  with 
a  41  x  81  grid  with  d  •  e  »  0.05  and  At  •  0.0125,  At  this  Reynolds 
number,  the  local  oscillations  in  the  comer  region  are  more  severe  than 
at  Re  «  9741. t>.  The  reversed  flow  regions  present  in  Figure  10  dissipate 
completely  by  t  ■  50.  Grid  convergence  studies  for  this  caso  indie* :e 
that  the  finite-difference  resolution  for  this  transformed  41  x  81  .rid 
is  inadequate  to  resolve  the  fine  scales  of  the  motion  present  in  tnc 
comer  region  at  early  time;  this  case  had  315  grid  points  in  the  comer 
region  dofined  by  0.9  <  r  <  1,  0  <  z  <  0.2.  Similar  calculations  for 
«  *  1,  Re  «  1  x  10$  developed  a  numerical  instability  at  t  ■  4.6;  the 
allowable  value  of  At,  given  by  liquation  (19),  approached  zero  due  to 
the  extreme  severity  of  the  local  oscillations  in  the  comer.  This  may 
indicate  the  development  of  a  physical  instability  at  this  high  Reynolds 
number. 

It  is  possible  that  the  reversed  flow  regions  observed  in  these 
calculations  are  related  to  a  physical  phenomenon  observed  in  Wcidraan's 
spin-up  experiments16  for  Re  ■  5.9  x  loV  He  conducted  experiments  for 
a  -  1.93  using  various  wall  acceleration  rates  and  found  that  for  wall 
accelerations  >  4  rad/s*  a  ’'turbulent  column"  formed  along  the  sidewall 
at  early  time  and  then  it  eventually  disappeared,...  "leaving  an  entirely 
laminar  approach  to  solid  body  rotation."  Although  the  largest  Reynolds 
number  used  in  our  calculations  is  much  lower  than  5.9  x  1(P,  one  can 
speculate  that  the  local  reversed  flow  regions  present  in  the  numerical 
calculations  at  moderate  Reynolds  number  are  manifestations  of  the 
observed  transient  "turbulent  column"  observed  by  Weidman.  Additional 
flow  visualization  experiments  are  required  to  investigate  this  further. 


VII.  EKMAN  LAYER  MASS  FLOW  AND  COMPATIBILITY 
CONDITION  DURING  SPIN-IJP 

An  accurate  description  of  the  Ekman  layer  radial  mass  flow  is  needed 
to  establish  an  appropriate  "compatibility  condition"  for  use  in  the 
Nedemcycr  spin-up  model The  Wedemeyer  model  has  been  used  by  many 
invest igators* * ‘ l6* 1 7» 1 8  to  study  spin-up  from  rest.  In  this  model 
the  flow  is  split  into  two  parts:  the  Ekman  layer  flow  and  the  remainder 

l6P.  D ,  Ueidman,  "On  the  Spin-Up  and  Spin-Doirt  of  a  Rotating  Fluid,  Part 
2.  Measurements  and  Stability,"  J.  Fluid  Meoh.,  Vol.  77,  Pt.  4,  1976, 
pp.  709-736. 

,7P.  D.  Ueidman,  "On  the  Spin-Up  and  Spin-Dcun  of  a  Rotating  Fluid ,  Part 
1.  Extending  the  Vedemsyer  Model,"  J .  Fluid  Me oh. ,  Vol.  77,  Pt.  4, 
1976,  pp.  686-70$. 

Goller  and  H .  Ranov,  "Unsteady  Rotating  Flam  in  a  Cylinder  vith  a 
Free  Surface,"  J .  Basic  Ena. .  Trans .  ASMS,  Vol.  90,  Series  D ,  December 
1968 ,  pp.  446-464. 
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of  the  flow,  culled  the  core  flow.  A  partial  differential  equation  for 
the  core  flow  is  derived  from  an  order  of  Magnitude  analysis  of  the 
Navier-Stokes  equations.  In  order  to  solve  this  equation,  Nedemeyer 
postulates  a  compatibility  condition,  or  functional  relationship  between 
the  radial  and  azimuthal  volocity  components  in  the  core,  to  approximate 
the  coupling  between  the  core  flow  and  the  flow  in  the  Ekman  layers. 

The  rodial  mass  flow  rate  in  the  EkMan  layer,  at  a  given  radial  position, 
prescribes  the  core  radial  velocity  by  conservation  of  mass.  Nedemeyer 
developed  his  linear  compatibility  condition  by  interpolating  between 
known  results  at  t  ■  0  and  it  can  be  illustrated  as  shown  in  Figure 
11.  He  assumed  that  this  condition  was  valid  as  long  as  Re  <  3  x  10s 
and  the  F.kman  layer  remained  laminar.  Several  invest igators9* 16 »17* 18 
have  used  Rogers  and  Lance's  numerical  solutions19,  for  the  laminar 
boundary -layer  flow  on  an  infinite  rotating  disk,  in  an  attempt  to  con¬ 
struct  a  more  accurate  "non-linear”  compatibility  condition.  These 
solutions,  for  various  ratios  of  the  outer  flow-to-disk  rotation  rate, 
have  been  utilized  to  produce  the  non-monotonic  compatibility  condition 
shown  in  Figure  11.  Neidaan16*17  used  the  Rogers  and  Lance  compatibility 
condition  together  with  the  Nedemeyer  model  and  found  that  the  non¬ 
monotonic  behavior  of  the  compatibility  condition  led  to  unrealistic 
double-valued  solutions  for  azimuthal  velocity.  Since  it  is  apparent 
from  Wcidman's  results  that  the  Rogers  and  Lance  compatibi lity  condition 
may  be  inappropriate  for  the  spin-up  problem,  we  have  attempted  to  use 
the  present  numerical  technique  to  quantify  the  Ekman  layer  mass  flow 
rate  during  spin-up  and  investigate  the  degree  of  applicability  of  both 
the  Medemeyer  and  Rogers  and  Lance  conditions. 

In  order  to  detenaine  the  instantaneous  outward  radial  mass  flow  at 
a  particular  radial  position  along  the  endwall  we  had  to  adopt  a  defini¬ 
tion  for  the  "boundary- layer"  edge  in  the  Navier-Stokes  calculations, 

Ne  define  this  edge,  6,  to  be  the  last  axial  position  away  from  the  end- 
wall  where  u  passes  through  zero.  The  numerical  results  show  u  passes 
through  zero  only  once  in  the  interior  at  a  particular  radial  position 
unless  there  is  a  temporary  reversed  flow  region(s)  present  at  that 
radial  position.  A  typical  radial  velocity  profile  obtained  from  the 
Navier-Stokes  solutions  is  shown  in  Figure  12  end  compared  with  a 
corresponding  Rogers  and  Lance  boundary  layer  result  for  a  non-rotating 
outer  flow  (the  Von  Karman  problem);  the  core  flow  above  the  Ekman  layer 
was  not  rotating  for  t  *  8.3  and  r  ■  0.76  in  the  present  calculations. 

The  two  results  in  Figure  12  are  almost  identical  as  they  should  be 
according  to  the  Nedemeyer  model.  The  small  differences  are  probably 
due  to  the  fact  that  u  does  not  asymptotically  approach  zero  at  the 
edge  of  the  Ekman  layer  for  spin-up  in  a  finite  cylinder.  The  Rogers 
and  Lance  calculations,  on  the  other  hand,  impose  this  asymptotic 
behavior  as  a  boundary  condition. 


l9W.  H.  Rogtre  and  G.  H.  Lanoe,  "The  Rotationally  Symmetric  Flai  of  a 
Vieooue  Fluid  in  the  Preeenoe  of  an  Infinite  Rotating  Disk,  *  J,  Fluid 
Heoh Vol.  ?t  Pt.  4 ,  April  I960 ,  pp.  617-931, 
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lor  purposes  of  the  present  comparison,  the  nondlmonsional  llkmun 
laser  radial  mass  flow  rate,  A,  is  determined  from  the  numerical  results 
us  inn  the  relation 


m  =  ¥  /Ke/r‘;  (21) 

O 

where  represents  the  value  of  t  at  t  »  i,  for  a  specified  radial 
position.  Hie  value  of  v  at  z  =  6  is  used  to  define  the  local  value  of 
v/r  for  comparison  with  the  coinpat ibi 1 ity  conditions  in  Figure  11.  Both 
't  <  and  are  obtained  from  the  numerical  results  by  linear  interpoint ion. 

Figure  13  shows  typical  numerical  results  for  m  as  a  function  of  v/r. 
Values  for  a  *  2,  He*  9741.0  are  plotted  for  three  radial  positions  and 
various  values  of  t,  indicated  by  the  numbers  adjacent  to  several  points. 
According  to  these  calculations  the  F.kmnn  layer  forms  during  approxi¬ 
mately  the  first  cylinder  rotation,  t  *  2n;  thereafter,  it  monotonical ly 
decays  as  t  and  v/r  increase.  The  data  in  Figure  13  were  obtained  using 
results  from  21  x  21  and  41  x  41  grids,  together  with  Richardson  quad¬ 
ratic  extrapolation*,  to  approximate  results  for  zero  grid  size.  The 
validity  of  this  approach  was  partially  affirmed  by  approximating  m  by 
the  first  two  terms  of  e  series 


m . 
t 


♦  h  (An.)* 


(22) 


where  h  and  i  are  unknown  constants  and  ft-  represents  the  value  of 
m  obtained  with  a  particular  grid  size  An,*.  Numerical  results  for  AB^  * 
An.*  «  0.100,  0.050  and  0.025  were  used  to  determine  the  three  constants 
in  l.quation  (22)  for  representative  r  and  t.  The  results  show  that 
1.0  <  i  <5.7,  indicating  that  the  actual  variation  is  no  more  extreme 
than  that  given  by  Richardson's  result  with  t  ■  2.  The  results  in 
Figure  13  arc  typical  of  those  for  calculations  performed  over  the 
parameter  range  1  <  a  <  4.4  and  405  <  Re  <  5  x  10*.  The  complete  set 

of  numerical  results,  together  with  grid  convergence  studies  performed 
for  i  =  1 ,  Rc  »  1000,  9741.0  and  5  x  10*.  supports  the  following  con¬ 
clusions: 


a.  the  F.kmnn  layer  mass  flow  at  a  given  radial  position  decreases 
monotonical ly  as  v/r  and  t  increase  for  t  >  2s; 

b.  there  is  no  unique  compatibility  condition  that  is  valid  in  an 
exact  sense  for  all  a,  Re,  r  and  t; 

c.  a  new  compatibility  condition  can  be  constructed  that  gives  a 
bettor  approximation  to  the  present  numerical  data  than  either  the 
Rogers  and  Lance  or  the  Wedcmcyer  compatibility  condition. 
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The  data  for  1  <  .»  «i  4,4  and  405  <  Re  <  5  x  I04  have  hocn  used  to 
develop  a  new  compatibility  condition  and  to  determine  its  error.  It 
was  constructed  by  expressing  the  average  value  of  i  lit  r  *  0.5  as  a 
function  of  (v/r).  A  simple  monotonic  function  that  approximates  the 
numerical  data  is 


( .443/9) ( 1ft  (v/r)3  -  24  (v/r)2  ♦  9J  for  0  <  (v/r)  <  0.75 
m  a  (25) 

(.445)  (l  *  (v/r)  |  for  0.75  <  (v/r)  <  1. 


It  is  plotted  in  figure  14  and  compared  with  the  Wetjemeyer  and  Rogers 
and  lance  conditions.  In  general,  the  calculated  m  is  larger  than  pre- 
uicted  by  liquation  (25)  for  r  <  0.5  and  smaller  for  r  >  0.5;  this  trend 
is  apparently  due  to  a  radial  variation  in  m  caused  by  the  presence  of 
the  sidewall.  The  curve  described  by  Equation  (25)  falls  between  the 
Wedemeyer  and  Rogers  and  luince  curves  for  0  <  v/r  <  0.75  and  is  coincident 
with  the  Wedemcyer  curve  for  0.75  <  v/r  <  1;  it  has  a  continuous  first 
derivative  over  the  interval  (0,1 J.  The  shaded  band  represents  the 
maximum  "scatter"  present  in  the  numerical  data  for  t  >  2s,  l  <  c»  <  4.4, 
405  <  Re  <  5xl04.  for  0  <  t  <  2s  the  mass  flow  is  less  than  that  in¬ 
dicated  by  the  shaded  band  for  a  given  (v/r);  this  early  time  behavior  is 
not  described  by  Equation  1*3 i  and  is  outside  the  scope  of  the  Ncdemeyer 
theory.  Note  that  Equation  (23)  is  not  fitted  to  the  center  of  the 
shaded  band;  rather,  it  is  purposely  fitted  to  the  r  »  0.5  data  to  weight 
the  results  toward  the  center  of  the  endwall  and  attempt  to  make  it 
equally  valid  for  r  <  ,5  and  r  >  .5.  A  curve  fitted  to  the  center  of 
the  shaded  band  would  weight  the  results  toward  the  sidewall  where  the 
radial  variation  in  4i  is  most  pronounced.  This  would  result  in  a  rela¬ 
tively  poor  approximation  to  m  near  the  center  of  the  cylinder.  It 
should  be  noted  that  the  adoption  of  the  compatibility  condition  given 
by  Equation  (23)  leads  to  errors  in  the  nondimcnsional  mass  flow  rate 
predicted  at  a  given  r  and  t  (and  hence  the  core  radial  velocity  in  the 
Wcdcmeyer  model)  of  as  much  as  approximately  20\  of  the  maximum  value, 
0.443.  This  rather  large  error  is  inherent  in  any  compatibility  con¬ 
dition  that  neglects  the  radial  variation  of  mass  flow.  We  expect  that 
the  compatibility  condition  given  by  Equation  (23)  should  be  equally 
valid  for  a  >  4.4  and  Re  >  5xl04  as  long  as  the  Ekman  layer  remains 
laminar.  A  compatibility  condition  appropriate  for  a  turbulent  Ekman 
layer  cannot  be  developed  using  result.*  from  the  present  numerical 
procedure. 
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VIII.  CONCLUSIONS 


A  prod l c t or - cor roc tor  milt iplc- iteration  mothod  has  been  combined 
with  the  Causs-Seidel  Iteration  technique  to  produce  an  accurate  and 
efficient  numerical  procedure  for  solving  the  unsteady  Navier-Stokcs 
equations.  Test  calculations  for  spin-up  in  a  cylinder  were  shown  to  be 
consistent  with  previous  calculations  and/or  experimental  measure¬ 
ments.  Computations  carried  out  for  0.5  <  a  <  4.4  and  205  <  Re  <  50,000 
showed  that  coordinate  transformations  could  be  used  to  simultaneously 
resolve  details  of  both  the  interior  and  boundary  layer  flows  using  a 
moderate  number  of  grid  points.  These  calculations  demonstrated  the 
presence  of  inertial  oscillations  and  temporary  reversed  flow  regions 
along  the  sidewall  during  spin-up  from  rest.  Computational  experiments 
have  been  performed  to  determine  the  radial  mass  flow  in  the  cndwall 
likman  layers  during  spin-up  and  investigate  the  accuracy  of  compatibility 
conditions  used  in  the  Wedemeyer  spin-up  model.  A  new  compatibility 
condition  has  been  developed  based  on  the  results  of  this  numerical 
study . 
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APPliNDIX  A 

TRANSFORM!.!)  EQUATIONS  AND  BOUNDARY  CONDITIONS 

The  Navier-Stokes  equations,  liquations  (1)  -  (5).  and  boundary 
conditions,  liquations  (11),  are  transformed  from  the  (r.t)  physical 
plane  to  the  (R,n)  computational  plane  using  the  analytical  transforma¬ 
tions  given  by  liquations  (12).  The  transformed  governing  equations  arc 

Yt  -  y8  1^"'  (*rr  *  Vr)  '  6r  nz  V'1 

♦  Yn  IRO*1  n„  ♦  8r  n2  fa/r) 

♦  Re*1  |  (B*)2  Ypa  ♦  (n*)2  .  (A*J 


‘•t  •  lRc  1  <»rr  4  *r/r>  *  8r  nz  Vr| 

4  :n  |Rc*!  n:z  *  Br  ni  Vr) 

y  <y  J 

♦  Re*  |  ( B ' )  “  ^  Cnni  ♦  2  YYn  \/r 

*  (.  l-Re'l/r2  ♦  n.  y  /r‘|  ,  iA*2) 

**  n 


where 


’  'W*  \n  "  rt/(Hr)2  *  V  V(V‘ 

♦  i l/(r  Br>  -  Brr/(er)‘l  ; 


(A- 3) 


B  *  2b/ ((b  ♦  r)  (b  -  r)  In  |  (b  ♦  P/ (b  -  1) ) )  ,  (A-4) 

r 


6  •  2rS  /((b  ♦  r)  (b  -  r)J, 

rr  r 


(A-S) 


«  2c/(a(c  ♦  z/a  -  l)(c  -  i/a  ♦  1)  In  [(c  ♦  l)/(c  -  1)).  (A-6) 


n  -  2(z/a  -  1)  n„/|a(c  ♦  i/a  -  1)(*T  -  r/a  ♦  1)), 
Z  * 


(A-7) 
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AmiNUlX  B 

(INI  Tit  -  D  2  FFl: RINCE  LQUAT I ONS 

The  finite-difference  approx i wit  ions  given  by  liquations  (15)  form 
the  basis  of  the  finite-difference  equations  representing  liquations  (A-l) 
and  (A-2).  Uncoupled  linear  difference  equations  arc  constructed  with 
the  I’CMI  method  by  approximating  T„  and  in  Equations  (A-l)  and  (A-2) 
and  >  and  >  in  Equation  (A-2)  as  known  constants,  specified  from  cither 
the  ext rnpolatcd  value  or  the  value  calculated  at  the  previous  iteration 
level.  The  coupling  between  equations  and  the  nonl ineari t ies  of 
individual  terms  are  approximated  by  the  multiple-itorution  process 
which  updates  T&,  tf ,  >  and  >  at  each  subsequent  iteration  level.  In 
the  present  difference  procedQre,  station  (i*l,j,k)  in  Figure  3  is  calcu¬ 
lated  from  known  information  at  stations  (i*l,j,k-i),  (i»l,j,k*l), 

(i.j.k)  and  (i-l,j.k).  Stations  (i*l,j-l,k)  and  (i*l,j*l,k)  are  treated 
as  unknowns,  resulting  in  a  tridiagonal  system  of  equations  for  liquations 
(A-l)  and  (A-2).  The  i -difference  equation  is  given  by  the  tridiagonal 
system 


m»  1 

alj  V»,j-1 


m+1 

'lj  Vl.j.k 


Clj  Vl.jM.k  *  Ulj 


(B-l) 


where 


ali  *  l(*rr  *  Vr)/Rc  ‘  flr  V,/r)/(2AB)  -  (Br/AB)2/Re.  (B  2) 


b,.  =  3/ (2At )  ♦  2  ( (£  /AS) 2  ♦  (n./An) 2 1/Rc, 
1  j  r 


c^  •  -  -  2  (Br/AB)  /Re, 


(B-3) 

(B-4) 


du  ■  <4 


(m) 


(") 


*  * \  Vr,/,2A'” 

•  •  *!:!., .k.,> 


(B-5) 
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The  r.-differencc  equation  is  given  by  the  tridiagonal  system 


2j  rUl,j-l.k  b2j  Ct*l.),k  2j  d2)  ’  (B  6) 


where 


«2j  "  *ij  ♦  ®r/(r  Ag  Re), 


(B-7) 


b2j  "  blj  4  <l/Rc  '  nt  1'n^r  ’ 


(B-8) 


c2j  ■  -  2  Br/(r  Re) , 


(B-9) 


d2j  ’  14  S.j.k  ‘  Vl.j,k)/l2At) 


♦  c(n)  -  Cl-) 

lt*l.j.k*l  S*l.j.k-1)  lnri/Re  ♦  n2  f-r  T6/r)/(2An) 

♦  Ulm!  ,  .  ♦  C(B!  .  .  (n  /An)2/Re 

» ♦  1 , .1 , k ♦  1  i ♦  1 ,  j  ,k- 1 )  K  i' 

♦  (vtn)  1  (y^  -  y^™^  )(n  /An)/r“' 

lVl.),k'  iYi*l.j,k*t  Yi*l.j,k-l,lVAn,'r  * 


(B- 10) 


The  finite-difference  approx i sat  ions  given  by  Equations  (14)  *011* 
the  basis  of  the  finite-difference  equations  representing  Equation  (A-3) . 
The  quantity  c  in  Equation  (A-3)  is  specified  from  either  the  extrapolated 
value  or  the  value  calculated  at  the  most  recent  t-iteration  level.  In 
the  Gauss-Scidcl  procedure  used  to  solve  Equation  (A-3),  station 
(i*l,j,k)  is  calculated  from  known  information  at  stations  (i*l,j-l,k), 
(i*l,j*l,k),  (i«l,j,k-l)  and  (i«l,j,k*l).  The  t-difference  equation  is 
given  by 


,(®*D  *  -(■♦!) 

l.j,k  *  j(Tiel.Jal,k  4  YiM,M,k)(VAe) 


4  (ti"!,J,kol  4  fiTl!j,k-l)(nx/AT*)4 


4  ^Tl.jM.k  ‘  ’!:;:j.l,k><‘rr  *  Vr>'<2  A6> 


(Cont'd) 
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(coin  Mi 


*  <’!?!. j.v.i  ■ 


(n) 


r  MM.j.kJ/,2(ftr/A‘<)‘  4  2lnz/An)‘|  » 


(B- 


where  t |  i  ^  and  ^  ^  arc  both  known  quantities  duo  to  the 

order  in  which  the  calculations  are  performed. 

The  difference  equations  shown  in  this  Appendix  are  solved  using 
the  numerical  procedure  described  in  Section  111. 
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1.1. ST  01-  SYMBOLS 


h,t,»o 

b.c.d.e 


u.v.w 


A6  ,  An 


cylinder  radius 

constants  in  liquation  (22) 

coordinate  transformation  constants 

cylinder  half>hoight 

Ikman  number  v/(Uc“)) 

time  subscript 

radial  subscript 

axial  subscript 

t- iteration  superscript 

nondimens ionul  Ikman  layer  radial  mass  flow  rate 
Y-iteration  superscript 

nondimcnsionnl  and  dimensional  radial  coordinate 

Reynolds  number  (•  U  a*/v] 

nondimensional  and  dimensional  time 

r,t»,z  nondimensional  velocity  components 

nondimensional  and  dimensional  axial  coordinate 

cylinder  aspect  ratio  (»  c/aj 

transformed  radial  coordinate 

nondimensional  and  dimensional  circulation 

llkman  layer  edge 

iteration  convergence  criteria 

nondimensional  and  dimensional  vorticity 

transformed  axial  coordinate 

azimuthal  coordinate 

liquid  kinematic  viscosity 

nondimensional  and  dimensional  stream  function 
final  cylinder  rotation  rate 
initial  cylinder  rotation  rate 

grid  sizes  in  6  and  n  coordinates 
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